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Computation of maximal local (un)stable manifold patches 
by the parameterization method 
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Abstract 

In this work we develop some automatic procedures for computing high order polynomial 
expansions of local (un)stable manifolds for equilibria of differential equations. Our method 
incorporates validated truncation error bounds, and maximizes the size of the image of the 
polynomial approximation relative to some specified constraints. More precisely we use that 
the manifold computations depend heavily on the scalings of the eigenvectors: indeed we 
study the precise effects of these scalings on the estimates which determine the validated 
error bounds. This relationship between the eigenvector scalings and the error estimates 
plays a central role in our automatic procedures. In order to illustrate the utility of these 
methods we present several applications, including visualization of invariant manifolds in the 
Lorenz and FitzHugh-Nagumo systems and an automatic continuation scheme for (un)stable 
manifolds in a suspension bridge problem. In the present work we treat explicitly the case 
where the eigenvalues satisfy a certain non-resonance condition. 
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algorithms, a posteriori analysis 
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1 Introduction 

Invariant sets are fundamental objects of study in dynamical systems theory. Sometimes we are 
interested in an invariant set which is a smooth manifold, and we seek a representation of a 
chart patch as the graph of a function or as the image of a chart map. Semi-numerical methods 
providing high order formal expansions of invariant manifolds have a long history in dynamical 
systems theory. We refer to the lecture notes of Sirno [I] , the historical remarks in Appendix B 
of the paper by Cabre, Fontich, and de la Llave [2], the manuscript of Haro [3], as well as the 
book by Meyer and Hall [0] for more complete discussion of this literature. 

The present work is concerned with algorithms for computing local stable/unstable mani¬ 
folds of equilibria solutions of differential equations, with validated error bounds. The methods 
employed here have some free computational parameters and we are especially interested in 
choosing these in an automatic way. We employ the parameterization method of mum in our 
computations. This method provides powerful functional analytic tools for studying invariant 
manifolds. The core of the parameterization method is an invariance equation which conjugates 
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a chart map for the local stable/unstable manifold to the linear dynamics given by the eigen¬ 
values (see for example ([3| in Section [2j . Expanding the invariance equation as a formal series 
and matching like powers leads to homological equations for the coefficients of the series. These 
homological equations are solved to any desired order, yielding a finite approximation. 

Given a finite approximate parameterization we would like to evaluate the associated trun¬ 
cation error. An important feature of the parameterization method is that there is a natural 
notion of a posteriori error, i.e. one can “plug” the approximate solution back into the invari¬ 
ance equation and measure the distance from zero in an appropriate norm. Further analysis 
is of course necessary in order to obtain validated error bounds, as small defects need not im¬ 
ply small truncation errors. When the invariance equation is formulated on a regular enough 
function space it is possible to apply a Newton-Kantorovich argument to get the desired bounds. 

A uniqueness result for the parameterization method states that the power series coefficients 
are unique up to the choice of the scalings of the (un)stable eigenvectors [5]. This freedom in the 
choice of scaling can be exploited in order to control the numerical properties of the scheme. For 
example by increasing or decreasing the length of the eigenvectors it is possible to manipulate 
the decay rates of the power series coefficients, and thus inffuence the numerical stability of the 
scheme. 

One of the main findings of the present work is that the bounds required in the Newton- 
Kantorovich argument (see the definition of the radii polynomials bounds in (20)) depend in 
an explicit way on the choice of the eigenvector scalings. This result leads to algorithms for 
optimizing the choice of eigenvectors scalings under some fixed constraints. The algorithms 
developed in the present work complement similar automatic schemes developed in [71 (for 
computer assisted study of periodic orbits) and are especially valuable in continuation arguments, 
where one wants to compute the invariant manifolds over a large range of parameter values in 
an automatic way. 


Remark 1.1. The optimization constraints referred to above can be chosen in different ways 
depending on ones goals. For example when the goal of the computation is visualization of the 
manifold it is desirable to choose scalings which maximize the “extent” of the manifold in phase 
space (i.e. maximize the surface measure of the patch). On the other hand when the eigenvalues 
have different magnitudes then it may be desirable to maximize the image of the manifold under 
the constraint that the ratios of the scalings of the eigenvectors are fixed (this is especially useful 
in “fast-slow” systems). In other situations one might want to optimize some other quantity 
all together. Whatever constraints one chooses, we always want to optimize while holding the 
error of the computation below some specified tolerance. The main point of the present work 
is that whatever the desired constraints, the explicit dependency of the bounds on the scaling 
facilitates the design of algorithms which respect the specified error tolerance. 


Remark 1.2. We fix the domain of our approximate parameterization to be the unit ball in C m 
(where m is the number of (un)stable eigenvalues, i.e. the dimension of the manifold) and vary 
the scalings of the eigenvectors in order to optimize with respect to the constraints. Another 
(theoretically equivalent approach) would be to fix the scalings of the eigenvectors and vary the 
size of the domain. However the scalings of the eigenvectors determine the decay rates of the 
power series coefficients, and working with analytic functions of fast decay seems to stabilize the 
problem numerically. 


Remark 1.3. In many previous applications of the parameterization method the free constants 
were selected by some “numerical experimentation.” See for example the introduction and 
discussion in Section 5 of [8], Remark 3.6 of [9], Remark 2.18 and 2.20 of m, the discussion of 
Example 5.2 in OH. Remark 2.4 of [12], and the discussion in Sections 4.2 and 6 of [T2]. This 
motivates the need for systematic procedures developed here. 
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Remark 1.4. The algorithms developed here facilitate the computation of local stable/unstable 
manifolds. Once the local computations have been optimized one could extend or “grow” larger 
patches of the local manifold using adaptive integration/continuation techniques. This is a 
topic of substantial research and we refer the interested reader to the survey article [[IB]. See 
also the works of nu na nsi izi usi nsi ed] and the references therein. Combining these 
integration/continuation algorithms with the methods of the present work could be an interesting 
topic for future research. 

Remark 1.5. In the present work we employ a functional analytic style of a-posteriori analysis 
in conjunction with the parameterization method of 0 Ei eg. Moreover the arguments are 
framed in classical weighted sequences spaces following the work of ED [22]. There are in 
the literature many other methods for obtaining rigorous computer assisted error bounds on 
numerical approximations of invariant manifolds. The interested reader should consult the 
works of liu HU EU E31EU ESI ESI E3 ESI 1221 EHl EH EU E31EU for other approaches and 
results. 

Remark 1.6. In recent years a number of authors have developed numerical methods based on 
the parameterization method in order to compute invariant manifolds of fixed and equilibrium 
points (e.g. see HU |35l [36] for more discussion). The parameterization method can also be used 
to compute stable/unstable manifolds associated with periodic orbits of differential equations 
pm GHEE], as well as stable/unstable manifolds associated with invariant circles/tori [39l HO). 
Indeed the parameterization method can be extended in order to compute the invariant tori 
themselves m , leading to a KAM theory “without action angle coordinates”. For more complete 
discussion of numerical methods based on the parameterization method we refer to the upcoming 
book [32]. For the moment we remark that the optimization algorithms developed in the present 
work could be adapted to these more general settings. 

Our paper is organized as follows. In Section [2] we present briefly the parameterization 
method and discuss its behaviour with respect to some specific changes of variable. In Section [3] 
we give a way to numerically compute an approximate parameterization and then address the 
issue of finding a rescaling that maximize the image of the parameterization, while verifying some 
a posteriori bounds that ensure (in some sense) the validity of the approximate parameterization. 
One possible way of proving the validity of the approximation is to use the ideas of rigorous 
computation, which we detail in Section [4] We conclude in Section [5] by presenting the results 
obtained with our method to compute maximal patches of local manifolds for several examples. 
The codes for all the examples can be found at ISU- 

2 The parameterization method 

In this section, we introduce the parameterization method for the stable manifold of an equilib¬ 
rium solution of a vector field. The unstable manifold is obtained by time reversal. 

2.1 Invariance equation for stable manifolds of equilibria of vector fields 

We consider an ordinary differential equation (ODE) of the form 


y' = g{y)i (i) 

where g : M n — > is analytic. Assume that p e M n is an equilibrium point, i.e. g(p) = 0, 

and assume that the dimension of the stable manifold at p is given by n s < n. Denote (A 14), 
1 < k < n s the stable eigenvalues (that is 3i(Afc) < 0, for k = 1 ,..., n s ) together with associated 
eigenvectors, and denote A = diag( Ai,..., A ns ). 
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We want to find an analytic parameterization of the local stable manifold at p. So we look 
for a power series representation 


m = E 

H>o 


a 


( dl \ 


i GM n % a a = 


\0nJ 

U n) / 


( 2 ) 


with the classical multi-indexes notations |a| = a\ + ... + a Us and 9 a = 9 ... 9n" s , and assume 
that the parameterization / conjugates the flow <p induced by g to the linear flow induced by A, 
that is 

f[e M e) =<p{t,m). 

Differentiating with respect to t and taking t = 0, we get that / satisfies the invariance equation 


Df{8)k9 = g{f{8)), 

and to get a well-posed problem we add the following constraints 

f(0)=p, D/(0) = (Vi ... V n ,). 


( 3 ) 


( 4 ) 


Endow C ns with norm ||$||c"<» = max{|0j,| : k = 1,... ,n s }, where | • | denotes the complex 
modulus, and using that norm, denote by B v C C ns the closed ball of radius v centered at 0. 
We look for a parameterization / which is analytic on a ball B v C C ris with v > 0. We call the 
image f[B v \ a patch of the local invariant manifold. 


Remark 2.1. If some of the eigenvalues happen to be complex-conjugate, say Ai = A 2 , ..., A 2 m-i 
A 2 )n, it is easier to consider a power series / with complex coefficients (i.e. with a a G C n ) and 
acting on 9 G C ns . We can then recover the real parameterization by considering, for 8 G M ns , 

freal(Q 1, ■ • ■ , 9 n J = /(#1 + z6*2, 8\ ~ • • • , #2m-l + l9 2 mi @2m— 1 ^2m.+l j • • • ; 9 ng ). 

See (43) for a more detailed explanation of this fact. To be general in the sequel of our presen¬ 
tation, we will assume that / is a complex power series. 


Remark 2.2. We say that there is a resonance of order a between the stable eigenvalues if 


«iAi + • • • + Otn s Aji — A j 


( 5 ) 


for some 1 < j < n s . If there is no resonance for any a G N ns then we say that the stable 
eigenvalues are non-resonant. Note that if |a| is large enough then a resonance is impossible. 

It is shown in [5] that if g is analytic then © has an analytic solution / as long as the 
eigenvalues are non-resonant. Moreover the power series coefficients of / are uniquely determined 
up to the choice of the scalings of the eigenvectors. This abstract result does not however 
provide explicit bounds on the size of the domain of analyticity B u for the parameterization: 
hence the need for a-posteriori validation of our numerical computations. We also note that 
if there is a resonance then the invariance equation can be modified so that we conjugate to 
a polynomial (instead of linear) dynamical system [Sj 01], and that the later work just cited 
implements computer assisted error bounds for the resonant case using the radii polynomial 
approach. Adapting the methods of the present work to the resonant case will be the topic of a 
future study. It is clear from the work of |5] that in the resonant case the Taylor coefficients of 
the parameterization are unique up to the choice of the eigenvector scalings. What remains to 
be checked is that in the resonant case the eigenvector scalings appear in the radii polynomials 
in an explicit way (as is the case in for non-resonant eigenvalues, see Section [4|. 
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2.2 Change of coordinates 

Assume that / is a power series of the form <§ satisfying <© and @ (therefore it is a local 
parameterization of the stable manifold at p). Now consider a change of coordinates in C Hs , 
defined by some invertible matrix T E M Ug (C), and the new power series 

m = f(re). 


Thanks to (|3]) , we have that 

Df(9)A6 = Df(F6)FA9 = g(f(T6)) = g{f{6)). ( 6 ) 

So if T is such that TA = AT, then / also satisfies the invariance equation ([3]), together with the 
slightly modified conditions 

f(o)=p, Df(o) = r (v± ... v na ). (7) 


Remark 2.3. From now on we assume that F = diag( 7 i,... , 7 n3 ), which is sufficient to have 
TA = Ar (it is also necessary if the \k are pairwise distinct). We also assume that the 7 \ 
are all real positive numbers and that coefficients 7 \ corresponding to two complex conjugates 
eigenvalues are equal. Taking 7 \ real is natural if all the eigenvalues are real (and / is therefore 
a real power series). On the other hand if there are some complex-conjugate eigenvalues, say 
Ai = A 2 , then the recovery of a real parameterization as explained in Remark 2.1 uses the 


fact that the corresponding eigenvectors V\ and V 2 also are complex-conjugate, and that this 
property is propagated to all the coefficients of the parameterization when recursively solving 
the invariance equation ([9|. By taking 71 and 72 real and equal, we ensure that this property 
is conserved after the rescaling (namely 71 Vj = 72 ^), so that we can still easily recover a real 
parameterization. Admittedly, we could relax this hypothesis and only assume that 71 and 72 
themselves are complex-conjugate, but we will not consider this possibility here. 


As announced, we now consider F = diag( 7 i,... , 7 ns ), where 7j E M+ for all 1 < i < n s , and 
/ defined as f(ff) = f(F9). The above discussion shows that / is a new parameterization of the 
local manifold, since it satisfies © and 0 . Besides, the Taylor expansion of / can be easily 
expressed in terms of the Taylor expansion of /. Indeed if we write / as 


f(0) = Y. 3a 0 a , 

M>o 


then the coefficients are given by 

a a = a Q 7 ", (8) 

where 7 = ( 71 ,... ,7 Ua ) and again standard multi-indexes notations. Therefore it is enough to 
find one parameterization / of the local manifold (or more precisely its coefficients a Q ) to get all 
the re-parameterizations / (at least those given by a diagonal matrix T) without further work. 
Let us introduce an operator acting on sequences to express this rescaling in a condensed way. 

Definition 1. Given 7 = ( 71 ,... , 7 n s ), we define C (acting on a) component-wise by 

£ a (a) = 7 Q a Q , V|a| > 0. 

Therefore, if a is the sequence of coefficients of the parameterization /, the sequence of 
coefficients of the parameterization / defined as above is given by £(a). 
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3 How to compute / and maximize the local manifold patch 


In this section we present a method to compute numerically a parameterization of the manifold 
(that is the coefficients (a a )) and then choose a proper rescaling 7 to maximize the corresponding 
image. We assume in the sequel that the nonlinearities in g are polynomials. Note that this not 
so restrictive as it might first seems, as techniques of automatic differentiation can be used in 
order to efficiently compute the (Taylor/Fourier/Chebyshev) series expansions of compositions 
with elementary functions. The authors first learned these techniques from Chapter 4.7 of [ 45] , 
but the interested reader should also refer to the discussion and references in [3 001 • 

Automatic differentiation is also a valuable tool for validated numerics, as polynomial non- 
linearities are often more convenient to work with than transcendental ones. Since elemen¬ 
tary functions of mathematical physics (powers, exponential, trigonometric functions, rational, 
Bessel, elliptic integrals, etc.) are themselves solutions of ODEs, these ODEs can be appended 
to the original problem of interest in order to obtain a new problem with only polynomial non- 
linearities (but with more variables and more equations). Moreover, in many computer assisted 
proofs it is the dimension of the underlying invariant object, and not the dimension of the em¬ 
bedding space, that informs the difficulty of the problem. We refer for example the book of m 
for a much more complete discussion of these matters. We also mention that automatic differ¬ 
entiation has been combined with the radii polynomial approach in |48| in order to compute 
periodic orbits of some celestial mechanics applications. 

Of course automatic differentiation is not the only method which can be used in order 
to replace a transcendental vector field with a polynomial one. Any method of polynomials 
approximation can be used. A detailed survey of the interpolation literature is far beyond 
the scope of the present work, however we mention the works of PI (50] where one can find 
implementation details and fuller discussion of the literature surrounding the use of Chebyshev 
polynomials to expand transcendental nonlinearities and obtain computer assisted error bounds. 
We also note that general purpose software exists for carrying out these kinds of manipulations, 
even with mathematically rigorous error bounds [511 521 [551 . 

3.1 Computation of the approximate parameterization 

Let / be a power series as in Q, assume g is a polynomial vector field of degree d given by 

g{y) = b py P ' b P G 

W<d 

and plug it into the invariance equation ©>• We obtain 

(«iAi + • • • + a ns X ns )a a 9 a = ^ 63 i ^ h d ( a/? ) Q > ( 9 ) 

|a|>0 \/3\<d \M>0 / H>0|/3|<(2 “ 

where again we use multi-indexes notations, a P = * ... * , and * denotes the 

Cauchy product. Notice that the two conditions in Q imply that the coefficients of order 0 and 
1 are the same on both sides of (|9]). There are several ways to obtain an approximation of the 
coefficients (a Q )| a |>2 so that (J9J) is satisfied, one of them being to compute them recursively for 
increasing |a|. Here we present another method, which fits naturally with the ideas of rigorous 
computations exposed later in the paper. We define the infinite dimensional vector a = (aa)| a |>o 


6 


and the operator F, acting on a component-wise by 


F a (a ) 


' ao~P, 

a ei ~ 

(aiAi+ 


T (%n s An s ) ®a 


E h M 

\P\<d 


if a = 0, 
if at = ei, V 1 < i < n s 

V lal > 2. 


Finding a solving ([9| and the additional conditions Q is equivalent to solve 


F(a) = {F Q (a)}| a |> 0 = 0. (10) 

Given 7 = (71 ,..., 7n s ), finding a rescaled parameterization (that is solving © and 0 ) can 
also be expressed as finding the zero of the function F. which is defined the same way as F 
except for the indices |a| = 1: 


F a {a) = < 


ao — p, if a = 0, 

a ei - 7 if a = ei, VI < i < n s 

(aiAi + ... + a Us A n J a a - ^ bp (a^) , V \a\ > 2. 


( 11 ) 


|<d 


Notice that the discussion in Section 2.2 shows that F(a) = 0 if and only if F(C(a)) = 0. 


Remark 3.1. Since ao and the a ei are fixed by the additional conditions Q, we could also 
consider them as parameters and define F as (F a )| a |> 2 , acting only on (a a ) | a |> 2 - We do this 
for the examples of Sections |5.l| and |5.2[ but we keep the above definition of F and F when we 
use rigorous computation (Section [4] and example in Section 5.31, because it allows for a simpler 
presentation. 

Now we fix an integer N and define the truncated operator = (F a )ii <iV , acting on a 
truncated sequence a [7V] = (aa)| Q |<jv, by 


F a {aW) = 


o-o ~ P, 

0>ei Vi; 

(«iAi F • • • F oi ns \ ns ) Q-a ^ ) bp > 


if a = 0, 
if a = ei, V 1 < % < n s 
V 2 < lal < N. 


■|<d 


Since the problem is now finite dimensional, we can use Newton’s method to compute an ap¬ 
proximate zero of fF] . the rest of this paper, a will denote such an approximate solution 
completed with 0 for |a| > N. See Section [ 5 ] for explicit examples. Also note that the only 
property that really matters concerning the approximate parameterization a is that a a = 0 for 
all |c*| > N. As long as it satisfies this property, everything in the sequel will work, even if 
a was obtained in a different fashion than the one we just presented (for instance by solving 
inductively a finite number of homological equations). 

Remark 3.2. Taking N larger leads to a better approximation but at the expense of compu¬ 
tational cost, so its choice depends on how precise an approximation you need, and how much 
computational resources you have. 


3.2 Maximizing the image of the parameterization 

Now that we have an approximate parameterization, we focus on maximizing the image of the 
corresponding manifold, while checking that our approximation is still valid. The power series 
/ given by © is now considered as 

C n 


7 










for some v > 0. One approach in getting the largest possible image of / would be to maximize 
the v for which © is valid on B u . We give in Definition [2] and Definition [3] two different 
definitions of parameterization validity. 


Remark 3.3. For reasons of numerical stability, we always consider the parameter space B u for 


v = 1 and instead use the 7 introduced in the reparameterization of Section 2.2 as a parameter. 
Indeed, assume that the parameterization / is valid on B Ul for some u \, then proving that it 
still is on B U2 for a different 1/2 is equivalent to prove that f(6) = f(T9) = /(7 i$i,..., 7 n 3 0n s ) is 
valid on B Vl , with 7fc = ^ for all k. So we can always keep v = 1 and rather try to maximize 


the 7 fc for which / is valid on B\. 

Based on the previous remark, from now on, and for the rest of the present paper, we always 
fix v = 1 , and therefore drop all references to this parameter. 


Remark 3.4. If the eigenvalues are real and not all equal to the same value, it may be useful to 
consider different scalings for each direction, that is to take T = diag( 7 i,... , 7 na ) with different 
7 ^ rather than T = diag( 7 ,... , 7 ). Indeed in this work we aim at maximizing the surface of 
the manifold patch, but for some specific problem (a fast-slow system for instance), you may 
rather want to enlarge the manifold in one precise direction, in which case you should definitely 
consider different 7 *, for each k. 


In this paper we will use two different criteria to say that our parameterization is valid on 
B u . The first one is a numerical a posteriori estimate and the second is a rigorous validation. In 
order to measure the validity of a parameterization, we need to compute the norm of a sequence 
a = {o«}| Q |>o with a a G C n . For this, let us introduce the space 


n 1 def 
'-v 



G C and ||u||^i = y~) 
M>o 


I u a \iy\ a \ < 00 


Given a = (aQ,)| Q |> 0 , with a a 
product space 


/ a^\ 

U"7 


G C n , denote aW 



Then, consider the 


X = 



lx 


def 

= max 

l<i<n 


iW 


1 1 


< 00 


Remark 3.5. It will be usefull to represent linear operators acting on elements of X with 
(infinite) matrix/vector notations. To prevent any future ambiguity, let us precise the ordering 
we use in this paper for those vectors and matrices. Given a G X, we represent it as the (infinite) 
vector (na )| Q ,|> 0 where the a a are ordered by growing |a|, and by lexicographical order within 
the coefficients with same |a|. For instance, if n s = 2, 

a = (ao,0) «o,i; ^2,0) <2,1,1, «o,2, • • -) T ■ 


Notice that each a a is himself a vector of M n . For an (infinite) matrix M = (M a ^), a | |^| >0 
representing a linear operator on X, we use the same order for the rows and columns. Notice 
that each coefficient M a a is in fact a n by n matrix whose coefficient on row i and column j 
(i i) 

will be denoted as M y a 7 , so that 

n 

(Ma )<?= y 

|/ 3 | >0 j=l 
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We now give the two announced criteria to measure the validity of a parameterization. 

Definition 2. Fix a defect threshold s max > 0, a tr unca tion dimension N and an approximate 
solution computed using the method of Section [3. l[ Denote a = a^ N 1. We say that 

m = < e ° ( 12 ) 

\a\<N 


is defect-valid on B u if 


||F(a)||.Y < £n 


Equivalently, we say that a is defect-valid on B v if (13) holds. Given 7 = ( 71 ,., 
say that the rescaled parameterization C(a) is defect-valid on B v if 


(13) 

, 7 n s ), we also 


\\F(jr(d))\\ x < Sr 


(14) 


Remember that g is assumed to be polynomial, and so F is also polynomial, say of degree d. 
Since a a = 0 for |a| > N, then F a (a) = 0 for all |a| > d{N — 1) + 1. Thus the quantity ||F(a) ||x 
in (13) is only a finite sum and can be computed explicitly. 

Assume now that we have computed all the F a (a) for |a| < d(N — 1) (which can be quite 
long because of the Cauchy products coming from the nonlinearities). When we then consider 
some 7 = ( 71 ,... , 7 ns ) and the rescaled parameterization C(a), we get (using the fact that the 
nonlinearities are polynomial and the definition of the Cauchy product) that for all |a| > 0, 


F a {C{d)) = 1 a F a {d). 


(15) 


This way, the evaluation of ||F(£(a))||x for any rescaling is computationally cheap and thus it 
is rather straightforward to find the 7 for which the re-paranreterization C(a) gives the largest 
image of the manifold, while being defect-valid. Let us be a little more precise about this. 
Depending on our goal we use two different approaches. 


Method 1: We look for eigenvector scalings which maximize the surface measure, subject to 
the restriction that the rescaled parameterization is defect-valid. Therefore we find numerically 
a mesh of the compact set 

{7 € M"' | \\F(C(d))\\ x = s max } 

and then approximately compute the surface area of the image for each point of the mesh. We 
refer to Sections |5.1|and|5.2|for explicit examples in dimension 2 . 


Method 2: We want to emphasize some specific directions when computing the manifold. 
Therefore we fix some weights uj\ ,..., c o ns and consider only rescalings of the form 

7 = 7 (t) = (Lui,...,fu; n J. 


We then look for the largest t such that the rescaled parameterization is defect-valid. By doing 
so we obtain a manifold that stretches more in the directions with the largest weights. We refer 
to Sections 5.1 and 5.2 for explicit examples in dimension 2 where we stretch the manifolds in 
the slow direction. 


Remark 3.6. When there is only one stable/unstable eigenvalue (or a single pair of complex 
conjugate eigenvalues) then Method 2 reduces to choosing the largest possible scaling for the 
eigenvector (or for the complex conjugate pair of eigenvectors) so that the rescaled parameteri¬ 
zation is defect-valid. 
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Now we would like to present a different definition of validity of a parameterization, inspired 
by the field of rigorous computing. For this, we briefly review the ideas of rigorous computation. 


The idea is to reformulate the problem F(a) = 0 given in (10) and to look for a fixed point of a 
Newton-like equation of the form 

T(a) =a- AF(a) 

where A is an approximate inverse of DF(a ), and a is a numerical approximation obtained by 
computing a finite dimensional projection of F (in our case we called it F^). Let us explain 
how we construct A. Remembering that 

F a (a) = (aqAi + ... + a na \ ns ) a a - ^ bp (</) , V |a| > 2 


m<d 


we consider the following approximation for DF(a ) 

(DF^(a) 


ft = 


A 


N 


A f 

Si. A 


iV+l 


V 


•7 


where for each k > N, ft k is a finite bloc diagonal matrix, each of its diagonal block being of 
size n and of the form (aqAi + ... + ot ns \ na ) I n , where |a| = k and I n is the n by n identity 


matrix. In other words (see Remark 3.5) 

4 ( a a)\a\=k = ((“ 1^1 + ■ ■ ■ + Ot ns \ ns ) a a ) | Q | =fe . 
We then define an approximate inverse A of DF{a) as 

(AW 0 \ 


A = 


An 

0 An+i 


\ 


(16) 


•/ 


where A^ n 1 is a numerical approximation of DF^l(o) 1 while the A^ j are the exact 

inverses. We then prove the existence of a zero of F by using a contraction argument yielding 
the existence of a fixed point of T. A precise theorem is stated below, but just before that we 
need (given 7 = ( 71 ,..., 7 nJ) to define a rescaled operator 


T = I-AF 


(17) 


that we can use in a similar fashion to prove the existence of a zero of F. Remembering that 

F(C(a)) = CF(a) 

we have 


and therefore we consider 


DF(C(a)) = CDF(a)C~ l 

ft = Cft C- 1 and A = FACT 1 
-1 


(18) 


as approximations for DF(C(a)) and ^ DF(C(a))j respectively. 
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The rigorous enclosure of a solution follows by verifying the hypothesis of the following 
Newton-Kantorovich type argument. Our method, often called the radii polynomial approach , 
was originally developed to study equilibria of PDEs j53j and was strongly influenced by the 
work of Yamamoto [55]. The differences between the radii polynomial approach and the stan¬ 
dard Newton-Kantorovich approach are mainly twofold. First, the map F under study is not 
required to map the Banach space X into itself. This is often the case when the map F comes 
from a differential equation and results in a loss of regularity of the function it maps. Second, 
the approach does not require controlling the exact inverse of the derivative, but rather only an 
approximate inverse. This can be advantageous as controlling exact inverses of infinite dimen¬ 
sional linear operator can be challenging. For more details on the radii polynomial approach 
for rigorous computations of stable and unstable manifolds of equilibria, we refer to [43]. Given 
r > 0, denote by B r (a ) C X = (ll) n the ball centered at a E X of radius r. 

Theorem 3.7. Let 7 = ( 7 i,..., 7 n J E R” s . Assume that the linear operator A in ([16]) is 
injective. For each i = l,...,n, assume the existence of bounds Y = (y^\ ... ,Y^j and 
Z(r) = |zW(r),...,lW(r)J such that 


(f(C(a))-C(a)] {l) <Y® and sup (DT (£(a) + b)c ) W < Z®(r). (19) 

V ' 01 b,c£B r ( 0) ' ' 0 \ 


TO 


i 1 


If there exists r > 0 such that 

pW(r) ^ y(0 + zW(r-) - r < 0 , for all i = 1 ,..., 


i 1 


n 


( 20 ) 


then T : B r (C(a)) —>• B r (C(a)) is a contraction. By the contraction mapping theorem, there 
exists a unique a* E B r (C(a )) C X such that F(a*) = 0. Moreover, ||a* — C(a)\\x < r. 


As we see in Section 0 the bounds (r),..., P^ (r) given in ( | 20 [ ) can be constructed as 
polynomials in r and are called the radii polynomials. 

The statement of Theorem 3.7 is now used to define our second definition of validity of a 
parameterization, which is of course more costly than the first one but provides rigorous bounds. 


Definition 3. Fix a proof threshold r max , a truncation dimension N and an approximate solu¬ 
tion o. Given a numerical zero a of F and 7 = ( 71 ,..., 7 „ s ), we say that the parameterization 


C{a) is proof-valid on B u if there exists r > 0 such that condition (20) holds for some r <r 


In the next section we explain how the bounds Y and Z can be constructed so that they 
depend explicitly on the scaling 7 . Then, as for Definition [2j you only need to do the costly 
computations once for a (that is for 7 = ( 1 ,..., 1 )) and then the new bounds (and thus the new 
radii polynomials ?W) can be computed easily for any rescaling. Therefore the process of finding 
the rescaling 7 which maximizes the image of a manifold given by a proof-valid parameterization 
is also rather straightforward. We give in Section 5.3 an example of application where we 
explicitly compute the bounds Y and Z. 


4 Explicit dependency of the radii polynomials in the scaling 7 

In this section we construct the bounds Y and Z satisfying ( jTT)] ) with an explicit dependency on 
the 7 whose action is given by (J 8 [) . 
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4.1 The bound Y 

Proposition 4.1. The bound Y = (T^,... ,Y^) defined component-wise by 


Y® = (£AF(a)) W 


, Vl<i<n, 

ll ~ ~ 


satisfies (19). 


Proof. By definition of T, 


and we have that 


which yields the formula for Y. 


f (£(d)) - C{d) = CAF[a) 
AF(£(d)) = CAF(a), 


□ 


Remark 4.2. As previously mentioned, F a (a) = 0 if |a| > d(N — 1) + 1, and since A is of the 
form 


A = 



\ 


Tail 

/ 


where Tail is a diagonal matrix (see ©), then Y can be computed as a finite sum. Moreover, 
the Y bound can be expensive to evaluate, since it requires computing the Cauchy products 
involved in F(a), the matrix D which is the numerical inverse of the full and possibly large 
matrix HF^(a), and the product AF(a). However, once AF(a) is computed, we only need 
to do the component-wise multiplication defined by L and the finite sum corresponding to the 
il norm to get the bound Y for any rescaling 7 . Therefore, recomputing the bound Y for a 
different rescaling is cheap. 


4.2 The bound Z 

For the clarity of the exposition, we now assume that the nonlinearity of g (and thus of F) are of 
degree 2. We insist that the method presented here still holds for nonlinearity of higher degree 
(see for instance Eli). It is also worth mentioning that in the context of computing equilibria 
of PDEs in |54l 1561 57] and periodic orbits of delay differential equations in |58| . the bounds 
of the radii polynomials have been derived for general polynomial problems. Here, we decided 
that staying fully general would only obscure the point with notations, hence our restriction to 
quadratic nonlinearities. 

To compute the Z bound, we split DT(C(a) + b)c as 

DT{C{d) + b)c= (/ — AA^c + A (DF{C{a) + b)A^ c 

= (l-AA^c + A ( DF(C{a )) - A t ) c + D 2 F(C(a)){b, c) 

and we are going to bound each term separately. 


4.2.1 The bound Zq 

We start this section with a result providing an explicit formula for the l\ operator norm of a 
matrix. 
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Lemma 4.3. Let g n ,n a ,N = n 


'N+n s -1 


) and B E M Bnna N { C). For all c E (£*)”, 


(Fc [7V1 


(0 


ll 


<£* 


(id) 

B 


3=1 


,(j) 


f 1 ’ 


where 


TS-Ud) 

K □' = max 

0<|/3|<7V \ 1/ 


1 


E 

0<|a|<7V 


B 


{id) 

a,/3 




, V 1 < i,i < n. 


( 21 ) 


Remark 4.4. The matrix/vector product should be understood according to Remark 3.5 with 
Qn,n s ,N simply being the length of (c Q )| a i <Ar seen as a vector of complex numbers. Lemma 4.3 is 
just the computation of the matrix norm associated to the weighted vector norm defined on £j,. 

Proposition 4.5. Let B = Iun^n+i) - AW(DFW(a)) and 

2 

-1 


B = 


C [n] B . 


( 22 ) 


Let the bound Zq = (Zq X \ ..., Zq 1> ) defined component-wise by 




zf = V 1 < i < n. 


'dd) 


3=1 


Then 


((/ - AA f ) 


(0 


a 


< Zq\ V 1 < i < n, 


for all c such that ||c||y < 1. 


Remark 4.6. This bound can also be quite costly, because of the matrix-matrix multiplication 
required to get B. But again, once B has been computed, we only need to do the multiplication 

by the diagonal matrices associated do and to get B and then to compute the 

(i j) 

quantities KT to get the new bound for any rescaling. 

Proof. We start by noticing that 

I-AA* = c(l-AA^C~ 1 . 

Then by definition of A 1 and A, — AA^ cj =0 for all |a| > N and we have 


((/ - Ai*) 


(0 


ll 


C [N] (i nN(N+l) ~ D (DF^(o))| (£ [iV] ) 1 C W N () 




and Lemma 4.3 yields the formula for Zq. 


□ 


4.2.2 The bound Z\ 

In this section we will need two additional results. The first one is a quantitative statement that 
£\, is a Banach algebra and allows us to bound the nonlinear terms. 

Definition 4. Let «,»E £} v . We denote by u* v the Cauchy product of u and v, namely 

(u*v) a = Y u a-pvp, V |a| > 0, 

0</3<a 

where fi < a means fii < ct* for all 1 < i < n s and (a — /3)i = a, — ffi for all 1 < i < n s . 
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Lemma 4.7. 


V u, v G il, ||u * v\\ e i < \\u\\ e i \\v\\ £ i . 

The second one bounds the action of the (infinite) diagonal part of A. 
Lemma 4.8. Let d E X = (f/) n , such that d a = 0 for all |a| < N. Then 


(Ad) 


(0 


< 


4 N min |5R(A;)| 

l<Z<n s 


d w 


ft 


V 1 < * < n. 


These two lemma allow us to get the Z\ bound. 

Proposition 4.9. The bound Z\ = [z^\ ... ,z[ n ^ defined component-wise by 


z[ k) = 


E 


l<i<n 


,(fc) 

4 


+ E 


l< 2 j<n 


L(fc) 


(£(a)) 


(*) 


p- 


N min |3i(Ai)| 

l<i<n s 


Vl<Kn, 


satisfies 


(0 


ft 


< V 1 <i<n, 


(i (DF(£(a))-At) c ) 

/or all c such that ||c||^ < 1. 

Remark 4.10. This bound is not costly, as we only need to get C(a) from a (a component-wise 
multiplication) and then to evaluate a finite sum to get the £j, norm of C{a). 

Proof. We first prove the bound without rescaling (that is for 7 = (1,...,1)). By definition of 
At, (J^DF(a) - At) c) = 0 for all | a\ < N. For |a| > N, remember that the general expression 
for F is (for quadratic linearity) 

F a (a) = (cciAi + ... + a ns \ ns ) a a - Y b P ( a ^) > V M > 2 - 

l/?l<2 

Then, again by definition of At, the (aqAi + ... + a„ s X Us ) term cancels out in ^ DF(a ) — At) c) 
and what is left is 


((llP(a)-At)c) =-( Y, b hd } + 

\1 <i<n 


E bij ( 


(*) * AA 


a' ' * c 


V |a| > N, (23) 




where /3* must be understood as the multi-index with 1 at index i and 0 elsewhere, and ffij as 


the multi-index with 1 at indexes i and j, and 0 elsewhere. We then use Lemma 4.7 to get 


(K®-At)c) w < £ 


i (k) 

4 


V l<i<n 


40 


el 


+ E 

l<2j<n 


,(fc) 


i» 


ft 


41) 


H 


We now use Lemma 4.8 which yields 

El<i<n 


(a (l>F( 6) - At) c 


(fc) 


< 


•SF) 


Ad 


ft 


+ E 




bi k) 

aW 


c (t) 


Pi,j 


/l 

*”1/ 




ft. 


tV min |3?(Af)| 

l<L<n s 


and the formula for Z\ follows (in the particular case when 7 = (1,..., 1)), since we assumed 
that ||c|| y < 1. Now we want to get the general bound. First notice that 


A [DF(C(d)) - At) c = LA (d.F( d) - A f ) CT x c. 


(24) 
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Then, going back to (23) and using that a* £ i c = £ 1 (£(a) * c), we get for all |a| > N that 


((flF(o) - £~ x cj 


, 1<2<72 


£ b, (c-'c)" ’+ V 6^ (^‘((As))®. <"»)), 

l<2j<n 


(25) 

Then, since we only need to consider the action of the diagonal part of A (that is for |a| > N) 
we can commute A and £ in (24). Finally, applying £ to (25) the £ and £~ l cancel out and 
using again Lemma 4.8 we get the announced formula for Z\. □ 


4.2.3 The bound Z 2 


To get the last bound we need a last lemma, which is a combination of Lemma 4.3 and Lemma 4.8 
and thus provides a bound on the full action of A. 

Lemma 4.11. For any d E (£]f) n and for all 1 < i < n, 


(. Ad ) (< > 




d(<) „ +EA 


(hi) 

o\ ’ z— < 

" iAi 


d k‘) 


ft 


Proposition 4.12. 77ie bound Z 2 = ( 
1 


= ( ..., Z^) defined component-wise by 


v (k) I J- 7 -Ak,k) 

lfVmm|^(A0| ^ 


E 


dk) 

Pi,j 




AM 
l^k 1 


di) 

7 &,j 


, V 1 < K n, 


where 

satisfies 


^[N] = £[N] a [N] ^[JV]^ 


-1 


(AD 2 F{£{d)){b lC ) 
for all b and c such that ||6||y < 1 and ||c||_y < 1. 


M 


H 


<zf. 


Remark 4.13. The only costly part in this bound is to get D (and the quantities K^' 1 ), but 
we already needed to compute D for the Y bound. 

Proof. Again we prove the bound without rescaling first (that is for 7 = (1,..., 1)). Since we 
assume that F is quadratic, we get that 

D 2 F(d)(b,c) = - £ bfab®*^, 

l<ij<n 


(26) 


with the same conventions as in Section 4.2.2 for the ffi.j . Therefore, using Lemma 4.7 and since 


| Y < 1 and ||c||jy < 1 , 


(D 2 F(a){b,i 


(fc) 


< £ 

1 <i,j<n 


r (fc) 

41 - 


Lemma 4.11 then yields the formula for Z 2 (in the particular case when 7 = (1,...,1)). To get 


the general formula, we can compute 

AD 2 F(£{a))(b,c) = £AD 2 F(a)(£~ 1 b, £~ 1 c) 

= £A£~ 1 D 2 F(d)(b,c), 

where we used (£ _ 1 6 ) * (£ _ 1 c) = £~ 1 ( b * c) in (26). The infinite part of £A£~ l (for |a| > N) 
is the same as the one of A since the infinite part of A is diagonal. The only difference is that 

(£AC -1 )^ = £^D = D, which yields the formula for Zj- □ 
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4.3 Radii polynomials 

Let us sum up the results of the previous sections. 

Proposition 4.14. Given 7 = (71,..., 7n s ), we consider F defined as in We also consider 

d an element of X = such that {d a ) a = 0 for all |or| > N (in practice a numerical 

approximate zero of F) and the operator T defined by @ and @. Then the bound Y 

defined in Proposition \f. t\ and the bound 

Z(r) = (Z 0 + Z 1 )r+Z 2 r 2 , 


where Z$, Z\ and Z 2 are defined in Propositions f.5, and 4-12 respectively, satisfy the 


hypothesis (19) of Theorem 3.1 


3.7 is a quadratic polynomial. If there 


Then, for each 1 < i < n, P® defined in Theorem 
exists r* > 0 such that ?W(r*) < 0 for all 1 < i < n, then there exis ts an interval I = (ro,?’i) 
such that X*®(r) < 0 for all 1 < i < n and for all r € I. By Theorem 


3.7 


we know that for all 

r £ I, within a ball of radius r centered in C(a) their exists a unique local parameterization of the 
manifold. Moreover, if one wants to make this fully rigorous, a final step consists of computing 
the bounds Y and Z with interval arithmetic and then check, still with interval arithmetic, that 
pW( r ) is negative. 

Finally, if the goal is to get a proof-valid parameterization while having the largest possible 
image, we process as follows. We start by computing the bounds (and the associated radii 
polynomials) without rescaling. Then if X is empty, or if r rnax < ro, we can rescale a to C(a) by 
some 7 and then compute the interval X associated to the rescaled polynomials P^ (of course 
one should choose 7 *, < 1 ) but this time the computation of the coefficients of the polynomials, 
namely Y, Zq, Z\ and Z 2 , are much faster thanks to the formulas of the previous sections. 
Conversely, if ro is small compared to r max , we can rescale a to C(a) by some 7 , this time with 
7 /c > 1 larger and larger, which will give a larger and larger manifold patch associated to the 
rescaled parameterization, until we reach the limit of ro = r max . We explain more in detail how 


we do this on an example in Section 5.3 


5 Examples 

5.1 Defect-valid parameterizations for the Lorenz system 

As a first example, we consider the Lorenz system, given by the vector field 

1 a {y ~ x ) \ 

g(x,y,z) = px - y - xz , 

\ xy-fiz ) 

with standard parameter values : o = 10, (3 = | and p = 28. In this case it is well known 
that the origin has a two dimensional stable manifold. We detail on this example the method 
presented in Sections [ 2 ] and [3] to automatically compute a maximal patch of the local stable 
manifold at p = 0 . 

We start by recalling that the stable eigenvalues are 


Ai = ( a T 1 + \/{cr - l) 2 + 4erp 


and A 2 = —fi, 


together with the stable eigenvectors 

Vi = 


Ai+ct 

1 


and V 2 = 
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As explained in Section [ 2 j we look for a parameteriztion of the local stable manifold in the form 
of a power series /, which should satisfy the invariance equation 


Df(6) n A 0 J e = g(m). 


(27) 


together with the condition conditions 

f(0)=p and Df( 0 ) = (Vi V 2 ) . 

Notice that in this case the two stable eigenvalues are real and therefore we can directly work 


with a real power series defined on [—1, l] 2 . Expanding / into a power series, (27) rewrites as 


^2 (aqAi + a 2 A 2 )a Q <9" = ^ 


M>2 


where 


M>2 


Clcv — 


/ CT (aa' > ~ ^ 

pa^ — al 2) — (a^ * 9 a , 

^ (a^ * — Pcffl j 

(a^\ 

a>u 

\a>a ) 

So we set a 0 ,o = V , opo = U, ao,i = V 2 and define F = (F a )| Q |> 2 , acting on a = (a Q )| Q |> 2 , by 


/ 


F a (a ) = (aiAi + a 2 A 2 )a Q — 


a I offi - ai 1 ') 

4 1} — CLa' 1 — 

(a^ 1 ) * — fici'a 1 

Our goal is now to find a numerical zero a and then the rescaling 7 = ( 71 , y 2 ) so that the 
parameterization / defined as 

m = 


\ 


\ 


V \a\ > 2. 


/ 


L C «{a)0 c 

| a |>0 


V 9 e [—1, l] 5 


gives us the maximal patch of manifold, while checking (according to Definition]^ that \\F(£(a))\\x < 
Emax > which will ensure that C(a) is a good approximate parameterization. 

First we fix an integer N and consider a truncated version of F, that is 

^ [n1 = to 2 < w<k , 

for which we can numerically compute a zero a with Newton’s method. Then we fix an E ma x 
and use Method 1 described in Section [ 3 ] First we compute F(a), which can be done explicitly 
because by construction a a = 0 for any |a| > N, so for 7 = 1, F( l \a) = 0 for any |a| > N 
and for i £ {2,3}, F^ 2 \a) = 0 for any |a| > 2 N — 1 (because of the quadratic terms). Then 
we find numerically the curve in the plane ( 71 , 7 : 2 ) that corresponds to ||F’(£(a)) ||x = £max■ In 
our case, we took a sample of values of 71 and for each we looked for the largest 7 2 for which 
||-F(jC(a))||x < £max (as explained in Section [ 3 ] this doesn’t require much computations since 
the coefficient of F(a) are already known). Finally we compute the surface of the corresponding 
patch of the manifold along this sample and find its maximum. The results are displayed 
in Figure [TJ along with the results of similar computations for the unstable manifolds of the 
nontrivial equilibria, or “eyes,” of the attractor. 

By way of contrast we consider another parameterization of the local stable manifold at 


p but focusing on the slow direction given by V^. Therefore we apply Method 2 described in 
Section |sj we define the ratio g = 


and only consider rescalings of the form 7 = ( 71 , gji ). 
Then we simply find numerically the largest 71 such that the rescaled parameterization C(a) is 
defect valid, and obtain the results displayed in Figure [2j 
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(a) The curve of ( 71 , 72 ) for which \\F(C(a)) ||.y = e mal (b) The corresponding values of the surface area (again 
for the local stable manifold of the origin. for the local stable manifold of the origin). 



(c) Lorenz System: local stable manifold of the origin (with the rescaling maximizing 
the surface area, i.e. 71 = 1.7 and 72 = 0.68) and local unstable manifolds of the eyes. 
The unstable manifolds have complex conjugate eigenvalues, so we simply maximize 
the length of the eigenvectors. 


Figure 1: For each manifold we take a defect constraint of e-max = 10 5 . The order of the 
parameterizations is N = 50 for the eyes and N = 30 for the stable local manifold of the origin. 

5.2 Defect-valid parameterizations for the FitzHugh-Nagumo equations 

We consider the vector field given by 

/ v \ 

g(u, v, w) = ^ (sv + w — q + u 3 — (1 + a)u 2 + cru) , 

V f (u - (w) ) 
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20 


Figure 2: Lorenz System: the figure illustrates the results of maximizing the lengths of the 
stable eigenvectors of the origin subject to the constraint that the slow eigenvector is g = | Ai|/| A 2 I 
times longer than the fast eigenvector and that the defect is less than Emax = 10” 5 . The order 
of this parameterization is N = 50. Note that the resulting patch covers more of the slow stable 
manifold than the patch shown in Figure [TJ however the surface area is smaller. 


where 


cr = —, s = 1.37, A = 1, q = 0.001, e = 0.15 and C = 5. 


There are trivial zeros given by v = 0, w = ^ and u solution of the cubic equation 

u 3 — (1 + cr)u 2 + ^(7 + u — q = 0. 

We want to compute the stable local manifold at one of them: 

/0.003374970076610\ 

p ~ 0 

^0.000674994015322 j 

With the selected values of the parameters we have two real stable eigenvalues at this point p: 

Ai ~ -0.662724919921474 and A 2 ~ -0.184083645070452, 
with associated eigenvectors 


-0.576099055982516\ 


—0.966141520359494\ 


Vi ~ I 0.381795200742850 and V 2 ~ 0.177850852721684 . 

, -0.722732524787547 J ^-0.186921472344981 J 


In this case we also want to compute a parameterization of the local stable manifold at p focusing 
more on the slow direction given by V 2 . Therefore we again apply Method 2 and obtain the 
results displayed in Figure [3} 
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Figure 3: FitzHugh-Nagumo System: the figure illustrates the results of maximizing the lengths 
of the stable eigenvectors of the origin subject to the constraint that the slow eigenvector is 
q = |Ai|/|A 2 1 times longer than the fast eigenvector and that the defect is less than 10 -5 . The 
order of this parameterization is IV = 30. The red star indicates the location of the equilibrium. 
The local manifold illustrated here is not the graph of any function over the stable eigenspace, 
i.e. the parameterization follows a fold in the manifold. Note that the triangulation in the figure 
is an artifact of the plotting procedure for the manifolds, and not part of the parameterization 
computation. If a finer mesh is desired we simply evaluate the polynomial approximation at 
more points. It is not necessary to re-compute the parameterization. 


5.3 Proof-valid parameterizations for the suspension bridge equation 

We consider the vector field 

( V 2 + V\V2 \ 

, \ V 3 

= 

\-(3v 3 -vi) 

which is obtained after a change of variable when one looks for travelling waves in the suspension 
bridge equation (e.g. see [591 W ) 


d 2 u 


d A u 

dx A 


(e u — 1) . 


We are going to rigorously compute the stable manifold at 0 (for a given /3 £ (0, 2)), which 
is two-dimensional. The stable eigenvalues are A and A, where 


^ — ~ 2 — /3 + *-\/2 + P , 

and associated eigenvectors are given by 


A 

A 2 

vW 


and V 2 = Vi. 


( 28 ) 
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We then define F = (-Fa)| a |> 0 , acting on a = (a a )| a |> 0 G (f 1 ) 4 , by 


F a (a) = < 


o 0 - 0, 

ago — hi, 

00,1 - Vi , 

(an A T d 2 A)dcK — 


(77 + (a (1) * a (2) ) a \ 


(3) 
a« 

(4) 

aA 


if a = 0 , 
if a = ( 1 , 0 ) 
if a = ( 0 , 1 ) 

V lal > 2. 


V Oq ^ - /daa ) ) 


This time since the eigenvalues are not real, we consider complex parameterization a, i.e. a a 6 C 4 
for all a. Then we compute a numerical zero a with the method described in Section |3.1| To 
rigorously prove the existence of a nearby solution a we follow the ideas exposed in Section [3] 
and consider an operator T of the form 

T : a a — AF(a). 

The following infinite matrix should be a good approximation of DF(a) (at least for N large 
enough) 

(DF^ n \cl) 0 \ 

An 

0 An+i 


A f = 


\ 




where for each k > N, A^ is a 4 {k + 1) by 4 (k + 1) bloc diagonal matrix defined as 

\ 

((k-l)X + X)h 

Ak = 


( fcA /4 0 

((k — 1)A + A )/ 4 

0 


V 

with I 4 the 4 by 4 identity matrix. Therefore we define 


kXFJ 


A = 


(D 0 

M n 

0 Mn+i 

V 


•7 


where D is a numerical approximation of DF^(a) 1 while the Aik = A, 1 are exact inverses. 
We are now ready to compute the bounds Y and Z defined in Section [3] in order to apply 


Theorem 3.7 an prove the existence of a true parameterization near a. In practice, we first 
compute the bounds without rescaling (that is for 7 = (1,1)) and denote them simply Y and Z, 
and then we find the largest rescaling for which the parameterization is still proof valid. 

5.3.1 Computation of the bounds Y and Z, and of the radii polynomials 

Concerning the bounds Y and Zq, there is nothing to add or to specify to what was said in 
Section |4j We set, for 1 < i < 4 


yW = (AF(a)) W 


i 1 ’ 
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and 




(*) _ 


= E^; 

3 =1 


(hi) 
B ’ 


(2 ?) 

where the K B are defined as in Section 


4.2.1 


For Z\ and Z 2 we can specify the bounds of 


Section [4j because we now work with a specific non linearity. We get 


z[ l) = 


1 + 

a« 

+ 

a< 2 > 






th 


AW)| 


z{ 2) = 


AW)| 


7 ( 3) = 1 

1 Awor 


zj 4) = 


1 +p 

EWE' 


and 


z! 2 l) = 2 max ( K {1,1) 


' 1 ^ 

D ’ |5i(A)|IV ) 


Z {2) = 2I\ b ’ 1 \ 


Z ( 2 ] = 2K§’ 1) , 


E 4) = 2 K$ X) . 


Now we can consider, for all 1 < i < 4, the radii polynomial defined by 

p«(r) = y« + - l)r + z£V 2 

and we can try and look for a positive r such that P^\r) < 0 for all 1 < i < 4. 

Remark 5.1. F should be very small if a is a good approximative zero of F. Z$ should also 

be very small because B = In(n+i) — D(DF^(a)) and D is a numerical inverse of (DF^I(a)). 

2 

Finally Z\ can be made very small by choosing N large enough. Therefore the radii polynomials 
are of the form 

pb)( r ) = e _ (1 _ rj^ r -f Z^r 2 

where e could be made arbitrarily small if we could get an arbitrarily good approximation a and 
7 1 could be made arbitrarily small if we could take with an arbitrarily large N (and if we could 
numerically compute inverses of matrices with sufficient accuracy). So up to having sufficient 
computational precision we should always be able to find a positive r such that P^\r) < 0. 


5.3.2 Results 

For this problem we are interested in proving (rigorously and with and error bound r smaller 
than r max ) the largest possible patch of the stable manifold, for values of j3 between 0.5 and 2. 
Since we already computed the bounds Y. Zq. Z\ and Z 2 without rescaling, we can now easily 
compute the radii polynomial P for any rescaling, and so we look by dichotomy for the largest 
7 such that the rescaled radii polynomial P has a positive root ro which is less or equal to r max . 
Notice that the eigenvalues are complex conjugated for this problem and that is why we only 
consider uniform rescaling (i.e. 71 = 72 ). 


When (3 goes to 2, the real part of A goes to 0 (remember (28)) so we expect it to be harder 


and harder to compute the manifold when f} goes to 2. Indeed we observe that the largest 7 for 
which we are able to do the proof becomes smaller and smaller when /3 goes to 2 (see Figure [4]). 
The computations were made with N = 30, v = 1 and r max = 10 -5 for the proof. 


Remark 5.2. Another interesting point here is that a closer look at the bound Z$ shows why 
it is better to take v = 1. T he matrix B is supposed to be approximatively 0, and we want the 

to be as small as possible, but their definition 
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Figure 4: Maximal value of 7 for which we can still do the proof with r < r max , for different 
values of (5. The manifold computations are completely automated. 


show that there is a risk of numerical errors if v is too small or too large, hence our choice of 
always considering v = 1 . 

To speed up the process of redoing the proof after a rescaling, we kept track of the 7 
dependency in the bound Y and Z, and constructed the rescaled bound Y and Z based on the 
original ones. However by doing things this way we introduce in the Zq bound the same kind of 
instability that comes with taking v 1 (see (22)). If the Zq bound becomes too big, we could 
deal with it (at the expense of speed), by recomputing all the bounds without using the fact 
that they came from a rescaling and thus eliminating this numerical instability issue. 
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